####################################
# Descriptive statistics
####################################

rm(list=ls())

#sink("~/Dropbox/Crime Chile/11_replication/01_descriptive.txt")

library(Hmisc)
library(stargazer)
library(foreign)
library(ggplot2)
library(rdrobust)
library(readstata13)

################
# Prepare data 
################

# reada data
d_original = read.dta13("~/Dropbox/Crime Chile/11_replication/entire_sample_chile_2022january.dta")   
d = read.dta13("~/Dropbox/Crime Chile/11_replication/local_crime_data_chile_2022january.dta")   
names(d)

###########################
# number of observations
###########################

nrow(d)
nrow(d_original)
nrow(d)/nrow(d_original)

######################
# Averages
######################

first_type = (mean(d$homicide_rate, na.rm = T) * sum(!is.na(d$homicide_rate))) + 
             (mean(d$rape_rate, na.rm = T) * sum(!is.na(d$rape_rate))) +
             (mean(d$assault_rate, na.rm = T) * sum(!is.na(d$assault_rate))) 

second_type = (mean(d$theft_rate, na.rm = T) * sum(!is.na(d$theft_rate))) +
              (mean(d$robbery_rate, na.rm = T) * sum(!is.na(d$robbery_rate))) +
              (mean(d$robbery_surprise_rate, na.rm = T) * sum(!is.na(d$robbery_surprise_rate))) 

total = (mean(d$homicide_rate, na.rm = T) * sum(!is.na(d$homicide_rate))) + 
        (mean(d$rape_rate, na.rm = T) * sum(!is.na(d$rape_rate))) +
        (mean(d$assault_rate, na.rm = T) * sum(!is.na(d$assault_rate))) +
        (mean(d$theft_rate, na.rm = T) * sum(!is.na(d$theft_rate))) +
        (mean(d$robbery_rate, na.rm = T) * sum(!is.na(d$robbery_rate))) +
        (mean(d$robbery_surprise_rate, na.rm = T) * sum(!is.na(d$robbery_surprise_rate))) 

round(first_type/total,2)
round(second_type/total,2)

round((mean(d$homicide_rate, na.rm = T) * sum(!is.na(d$homicide_rate)))/total,3)
round((mean(d$rape_rate, na.rm = T) * sum(!is.na(d$rape_rate)))/total,3)
round((mean(d$assault_rate, na.rm = T) * sum(!is.na(d$assault_rate)))/total,2)
round((mean(d$theft_rate, na.rm = T) * sum(!is.na(d$theft_rate)))/total,2)
round((mean(d$robbery_rate, na.rm = T) * sum(!is.na(d$robbery_rate)))/total,2)
round((mean(d$robbery_surprise_rate, na.rm = T) * sum(!is.na(d$robbery_surprise_rate)))/total,2) 

###############################################################
# Subsetting
################################################################

h1 = rdbwselect(d$homicide_rate,d$margin,cluster = d$cluster)$bws[1]
h2 = rdbwselect(d$rape_rate,d$margin,cluster = d$cluster)$bws[1]
h3 = rdbwselect(d$assault_rate,d$margin,cluster = d$cluster)$bws[1]
h4 = rdbwselect(d$theft_rate,d$margin,cluster = d$cluster)$bws[1]
h5 = rdbwselect(d$robbery_rate,d$margin,cluster = d$cluster)$bws[1]
h6 = rdbwselect(d$robbery_surprise_rate,d$margin,cluster = d$cluster)$bws[1]

d1 = d[abs(d$margin)<=h1,]
d2 = d[abs(d$margin)<=h2,]
d3 = d[abs(d$margin)<=h3,]
d4 = d[abs(d$margin)<=h4,]
d5 = d[abs(d$margin)<=h5,]
d6 = d[abs(d$margin)<=h6,]

###########################
# Eligible sample
###########################

crime = data.frame(d$homicide_rate,
                   d$rape_rate,
                   d$assault_rate,
                   d$theft_rate ,
                   d$robbery_rate,
                   d$robbery_surprise_rate)

colnames(crime) = c("Homicide",
                    "Rape",
                    "Assault",
                    "Theft",
                    "Robbery",
                    "Robbery by surprise")

colnames(crime)

stargazer(crime,
          type = "text",
          colnames = FALSE,
          min.max = FALSE,
          title="", 
          summary.stat = c("mean", "sd","n"),
          digits=0, 
          rownames=FALSE, 
          align=TRUE, 
          float = TRUE, 
          float.env = "table", 
          table.placement = "H",
          header = FALSE, 
          no.space = TRUE,
          out="~/Dropbox/Crime Chile/11_replication/table1.html")

###########################
# Covariate comparison
###########################

capital_d_original = mean(d_original$capitalprovince, na.rm = T)
province_d_original = mean(d_original$codeprovince, na.rm = T)
region_d_original = mean(d_original$coderegion, na.rm = T)
area_d_original = mean(d_original$area, na.rm = T)
population_d_original = mean(d_original$logpop02, na.rm = T)
urban_d_original = mean(d_original$urbano02_p, na.rm = T)
literacy_rate_d_original = mean(d_original$alfab03, na.rm = T)
income_d_original = mean(d_original$ingresopercapita03, na.rm = T)
hdi_d_original = mean(d_original$idh03, na.rm = T)
hdr_d_original = mean(d_original$rankidh03, na.rm = T)
votes_d_original = mean(d_original$logpresvotes1999_1_all, na.rm = T)
right_d_original = mean(d_original$lavin99_p, na.rm = T)
left_d_original = mean(d_original$lagos99_p, na.rm = T)
invalid_d_original = mean(d_original$nulo99_p, na.rm = T)
blank_d_original = mean(d_original$blanco99_p, na.rm = T)
unemployment_d_original = mean(d_original$unemployment_2003, na.rm = T)
wage_d_original = mean(d_original$wage_2003, na.rm = T)
inequality_d_original = mean(d_original$inequality_2003, na.rm = T)

capital_d = mean(d$capitalprovince, na.rm = T)
province_d = mean(d$codeprovince, na.rm = T)
region_d = mean(d$coderegion, na.rm = T)
area_d = mean(d$area, na.rm = T)
population_d = mean(d$logpop02, na.rm = T)
urban_d = mean(d$urbano02_p, na.rm = T)
literacy_rate_d = mean(d$alfab03, na.rm = T)
income_d = mean(d$ingresopercapita03, na.rm = T)
hdi_d = mean(d$idh03, na.rm = T)
hdr_d = mean(d$rankidh03, na.rm = T)
votes_d = mean(d$logpresvotes1999_1_all, na.rm = T)
right_d = mean(d$lavin99_p, na.rm = T)
left_d = mean(d$lagos99_p, na.rm = T)
invalid_d = mean(d$nulo99_p, na.rm = T)
blank_d = mean(d$blanco99_p, na.rm = T)
unemployment_d = mean(d$unemployment_2003, na.rm = T)
wage_d = mean(d$wage_2003, na.rm = T)
inequality_d = mean(d$inequality_2003, na.rm = T)

capital_d1 = mean(d1$capitalprovince, na.rm = T)
province_d1 = mean(d1$codeprovince, na.rm = T)
region_d1 = mean(d1$coderegion, na.rm = T)
area_d1 = mean(d1$area, na.rm = T)
population_d1 = mean(d1$logpop02, na.rm = T)
urban_d1 = mean(d1$urbano02_p, na.rm = T)
literacy_rate_d1 = mean(d1$alfab03, na.rm = T)
income_d1 = mean(d1$ingresopercapita03, na.rm = T)
hdi_d1 = mean(d1$idh03, na.rm = T)
hdr_d1 = mean(d1$rankidh03, na.rm = T)
votes_d1 = mean(d1$logpresvotes1999_1_all, na.rm = T)
right_d1 = mean(d1$lavin99_p, na.rm = T)
left_d1 = mean(d1$lagos99_p, na.rm = T)
invalid_d1 = mean(d1$nulo99_p, na.rm = T)
blank_d1 = mean(d1$blanco99_p, na.rm = T)
unemployment_d1 = mean(d1$unemployment_2003, na.rm = T)
wage_d1 = mean(d1$wage_2003, na.rm = T)
inequality_d1 = mean(d1$inequality_2003, na.rm = T)

capital_d2 = mean(d2$capitalprovince, na.rm = T)
province_d2 = mean(d2$codeprovince, na.rm = T)
region_d2 = mean(d2$coderegion, na.rm = T)
area_d2 = mean(d2$area, na.rm = T)
population_d2 = mean(d2$logpop02, na.rm = T)
urban_d2 = mean(d2$urbano02_p, na.rm = T)
literacy_rate_d2 = mean(d2$alfab03, na.rm = T)
income_d2 = mean(d2$ingresopercapita03, na.rm = T)
hdi_d2 = mean(d2$idh03, na.rm = T)
hdr_d2 = mean(d2$rankidh03, na.rm = T)
votes_d2 = mean(d2$logpresvotes1999_1_all, na.rm = T)
right_d2 = mean(d2$lavin99_p, na.rm = T)
left_d2 = mean(d2$lagos99_p, na.rm = T)
invalid_d2 = mean(d2$nulo99_p, na.rm = T)
blank_d2 = mean(d2$blanco99_p, na.rm = T)
unemployment_d2 = mean(d2$unemployment_2003, na.rm = T)
wage_d2 = mean(d2$wage_2003, na.rm = T)
inequality_d2 = mean(d2$inequality_2003, na.rm = T)

capital_d3 = mean(d3$capitalprovince, na.rm = T)
province_d3 = mean(d3$codeprovince, na.rm = T)
region_d3 = mean(d3$coderegion, na.rm = T)
area_d3 = mean(d3$area, na.rm = T)
population_d3 = mean(d3$logpop02, na.rm = T)
urban_d3 = mean(d3$urbano02_p, na.rm = T)
literacy_rate_d3 = mean(d3$alfab03, na.rm = T)
income_d3 = mean(d3$ingresopercapita03, na.rm = T)
hdi_d3 = mean(d3$idh03, na.rm = T)
hdr_d3 = mean(d3$rankidh03, na.rm = T)
votes_d3 = mean(d3$logpresvotes1999_1_all, na.rm = T)
right_d3 = mean(d3$lavin99_p, na.rm = T)
left_d3 = mean(d3$lagos99_p, na.rm = T)
invalid_d3 = mean(d3$nulo99_p, na.rm = T)
blank_d3 = mean(d3$blanco99_p, na.rm = T)
unemployment_d3 = mean(d3$unemployment_2003, na.rm = T)
wage_d3 = mean(d3$wage_2003, na.rm = T)
inequality_d3 = mean(d3$inequality_2003, na.rm = T)

capital_d4 = mean(d4$capitalprovince, na.rm = T)
province_d4 = mean(d4$codeprovince, na.rm = T)
region_d4 = mean(d4$coderegion, na.rm = T)
area_d4 = mean(d4$area, na.rm = T)
population_d4 = mean(d4$logpop02, na.rm = T)
urban_d4 = mean(d4$urbano02_p, na.rm = T)
literacy_rate_d4 = mean(d4$alfab03, na.rm = T)
income_d4 = mean(d4$ingresopercapita03, na.rm = T)
hdi_d4 = mean(d4$idh03, na.rm = T)
hdr_d4 = mean(d4$rankidh03, na.rm = T)
votes_d4 = mean(d4$logpresvotes1999_1_all, na.rm = T)
right_d4 = mean(d4$lavin99_p, na.rm = T)
left_d4 = mean(d4$lagos99_p, na.rm = T)
invalid_d4 = mean(d4$nulo99_p, na.rm = T)
blank_d4 = mean(d4$blanco99_p, na.rm = T)
unemployment_d4 = mean(d4$unemployment_2003, na.rm = T)
wage_d4 = mean(d4$wage_2003, na.rm = T)
inequality_d4 = mean(d4$inequality_2003, na.rm = T)

capital_d5 = mean(d5$capitalprovince, na.rm = T)
province_d5 = mean(d5$codeprovince, na.rm = T)
region_d5 = mean(d5$coderegion, na.rm = T)
area_d5 = mean(d5$area, na.rm = T)
population_d5 = mean(d5$logpop02, na.rm = T)
urban_d5 = mean(d5$urbano02_p, na.rm = T)
literacy_rate_d5 = mean(d5$alfab03, na.rm = T)
income_d5 = mean(d5$ingresopercapita03, na.rm = T)
hdi_d5 = mean(d5$idh03, na.rm = T)
hdr_d5 = mean(d5$rankidh03, na.rm = T)
votes_d5 = mean(d5$logpresvotes1999_1_all, na.rm = T)
right_d5 = mean(d5$lavin99_p, na.rm = T)
left_d5 = mean(d5$lagos99_p, na.rm = T)
invalid_d5 = mean(d5$nulo99_p, na.rm = T)
blank_d5 = mean(d5$blanco99_p, na.rm = T)
unemployment_d5 = mean(d1$unemployment_2003, na.rm = T)
wage_d5 = mean(d5$wage_2003, na.rm = T)
inequality_d5 = mean(d5$inequality_2003, na.rm = T)

capital_d6 = mean(d6$capitalprovince, na.rm = T)
province_d6 = mean(d6$codeprovince, na.rm = T)
region_d6 = mean(d6$coderegion, na.rm = T)
area_d6 = mean(d6$area, na.rm = T)
population_d6 = mean(d6$logpop02, na.rm = T)
urban_d6 = mean(d6$urbano02_p, na.rm = T)
literacy_rate_d6 = mean(d6$alfab03, na.rm = T)
income_d6 = mean(d6$ingresopercapita03, na.rm = T)
hdi_d6 = mean(d6$idh03, na.rm = T)
hdr_d6 = mean(d6$rankidh03, na.rm = T)
votes_d6 = mean(d6$logpresvotes1999_1_all, na.rm = T)
right_d6 = mean(d6$lavin99_p, na.rm = T)
left_d6 = mean(d6$lagos99_p, na.rm = T)
invalid_d6 = mean(d6$nulo99_p, na.rm = T)
blank_d6 = mean(d6$blanco99_p, na.rm = T)
unemployment_d6 = mean(d6$unemployment_2003, na.rm = T)
wage_d6 = mean(d6$wage_2003, na.rm = T)
inequality_d6 = mean(d6$inequality_2003, na.rm = T)

mean_d_original = round(c(capital_d_original,province_d_original,region_d_original,area_d_original,population_d_original,urban_d_original,literacy_rate_d_original,income_d_original,hdi_d_original,hdr_d_original,votes_d_original,right_d_original,left_d_original,invalid_d_original,blank_d_original,unemployment_d_original,wage_d_original,inequality_d_original),2)
mean_d =  round(c(capital_d, province_d, region_d, area_d, population_d, urban_d, literacy_rate_d, income_d, hdi_d, hdr_d, votes_d, right_d, left_d, invalid_d, blank_d,unemployment_d,wage_d,inequality_d),2)
mean_d1 = round(c(capital_d1,province_d1,region_d1,area_d1,population_d1,urban_d1,literacy_rate_d1,income_d1,hdi_d1,hdr_d1,votes_d1,right_d1,left_d1,invalid_d1,blank_d1,unemployment_d1,wage_d1,inequality_d1),2)
mean_d2 = round(c(capital_d2,province_d2,region_d2,area_d2,population_d2,urban_d2,literacy_rate_d2,income_d2,hdi_d2,hdr_d2,votes_d2,right_d2,left_d2,invalid_d2,blank_d2,unemployment_d2,wage_d2,inequality_d2),2)
mean_d3 = round(c(capital_d3,province_d3,region_d3,area_d3,population_d3,urban_d3,literacy_rate_d3,income_d3,hdi_d3,hdr_d3,votes_d3,right_d3,left_d3,invalid_d3,blank_d3,unemployment_d3,wage_d3,inequality_d3),2)
mean_d4 = round(c(capital_d4,province_d4,region_d4,area_d4,population_d4,urban_d4,literacy_rate_d4,income_d4,hdi_d4,hdr_d4,votes_d4,right_d4,left_d4,invalid_d4,blank_d4,unemployment_d4,wage_d4,inequality_d4),2)
mean_d5 = round(c(capital_d5,province_d5,region_d5,area_d5,population_d5,urban_d5,literacy_rate_d5,income_d5,hdi_d5,hdr_d5,votes_d5,right_d5,left_d5,invalid_d5,blank_d5,unemployment_d5,wage_d5,inequality_d5),2)
mean_d6 = round(c(capital_d6,province_d6,region_d6,area_d6,population_d6,urban_d6,literacy_rate_d6,income_d6,hdi_d6,hdr_d6,votes_d6,right_d6,left_d6,invalid_d6,blank_d6,unemployment_d6,wage_d6,inequality_d6),2)

colnames = c("Capital",
               "Province",	  
               "Region",     
               "Area",      
               "Population",
               "Urban",   
               "Literacy",	
               "Income",     
               "HDI",    
               "HDR", 
               "Votes",	
               "Right-wing",   
               "Left-wing",   
               "Invalid",    
               "Blank",
               "Unemployment",
               "Wage",
               "Inequality")

############################################
# Standardized differences d original and d
############################################

sd1_1 = abs((mean(d_original$capitalprovince, na.rm = T) - mean(d$capitalprovince, na.rm = T))/ sqrt( (sd(d_original$capitalprovince, na.rm = T)^2 - sd(d$capitalprovince, na.rm = T)^2/2) ))
sd2_1 = abs((mean(d_original$codeprovince, na.rm = T) - mean(d$codeprovince, na.rm = T))/ sqrt( (sd(d_original$codeprovince, na.rm = T)^2 - sd(d$codeprovince, na.rm = T)^2/2) ))
sd3_1 = abs((mean(d_original$coderegion, na.rm = T) - mean(d$coderegion, na.rm = T))/ sqrt( (sd(d_original$coderegion, na.rm = T)^2 - sd(d$coderegion, na.rm = T)^2/2) ))
sd4_1 = abs((mean(d_original$area, na.rm = T) - mean(d$area, na.rm = T))/ sqrt( (sd(d_original$area, na.rm = T)^2 - sd(d$area, na.rm = T)^2/2) ))
sd5_1 = abs((mean(d_original$logpop02, na.rm = T) - mean(d$logpop02, na.rm = T))/ sqrt( (sd(d_original$logpop02, na.rm = T)^2 - sd(d$logpop02, na.rm = T)^2/2) ))
sd6_1 = abs((mean(d_original$urbano02_p, na.rm = T) - mean(d$urbano02_p, na.rm = T))/ sqrt( (sd(d_original$urbano02_p, na.rm = T)^2 - sd(d$urbano02_p, na.rm = T)^2/2) ))
sd7_1 = abs((mean(d_original$alfab03, na.rm = T) - mean(d$alfab03, na.rm = T))/ sqrt( (sd(d_original$alfab03, na.rm = T)^2 - sd(d$alfab03, na.rm = T)^2/2) ))
sd8_1 = abs((mean(d_original$ingresopercapita03, na.rm = T) - mean(d$ingresopercapita03, na.rm = T))/ sqrt( (sd(d_original$ingresopercapita03, na.rm = T)^2 - sd(d$ingresopercapita03, na.rm = T)^2/2) ))
sd9_1 = abs((mean(d_original$idh03, na.rm = T) - mean(d$idh03, na.rm = T))/ sqrt( (sd(d_original$idh03, na.rm = T)^2 - sd(d$idh03, na.rm = T)^2/2) ))
sd10_1 = abs((mean(d_original$rankidh03, na.rm = T) - mean(d$rankidh03, na.rm = T))/ sqrt( (sd(d_original$rankidh03, na.rm = T)^2 - sd(d$rankidh03, na.rm = T)^2/2) ))
sd11_1 = abs((mean(d_original$logpresvotes1999_1_all, na.rm = T) - mean(d$logpresvotes1999_1_all, na.rm = T))/ sqrt( (sd(d_original$logpresvotes1999_1_all, na.rm = T)^2 - sd(d$logpresvotes1999_1_all, na.rm = T)^2/2) ))
sd12_1 = abs((mean(d_original$lavin99_p, na.rm = T) - mean(d$lavin99_p, na.rm = T))/ sqrt( (sd(d_original$lavin99_p, na.rm = T)^2 - sd(d$lavin99_p, na.rm = T)^2/2) ))
sd13_1 = abs((mean(d_original$lagos99_p, na.rm = T) - mean(d$lagos99_p, na.rm = T))/ sqrt( (sd(d_original$lagos99_p, na.rm = T)^2 - sd(d$lagos99_p, na.rm = T)^2/2) ))
sd14_1 = abs((mean(d_original$nulo99_p, na.rm = T) - mean(d$nulo99_p, na.rm = T))/ sqrt( (sd(d_original$nulo99_p, na.rm = T)^2 - sd(d$nulo99_p, na.rm = T)^2/2) ))
sd15_1 = abs((mean(d_original$blanco99_p, na.rm = T) - mean(d$blanco99_p, na.rm = T))/ sqrt( (sd(d_original$blanco99_p, na.rm = T)^2 - sd(d$blanco99_p, na.rm = T)^2/2) ))
sd16_1 = abs((mean(d_original$unemployment_2003, na.rm = T) - mean(d$unemployment_2003, na.rm = T))/ sqrt( (sd(d_original$unemployment_2003, na.rm = T)^2 - sd(d$unemployment_2003, na.rm = T)^2/2) ))
sd17_1 = abs((mean(d_original$wage_2003, na.rm = T) - mean(d$wage_2003, na.rm = T))/ sqrt( (sd(d_original$wage_2003, na.rm = T)^2 - sd(d$wage_2003, na.rm = T)^2/2) ))
sd18_1 = abs((mean(d_original$inequality_2003, na.rm = T) - mean(d$inequality_2003, na.rm = T))/ sqrt( (sd(d_original$inequality_2003, na.rm = T)^2 - sd(d$inequality_2003, na.rm = T)^2/2) ))

sd_1 = round(c(sd1_1,sd2_1,sd3_1,sd4_1,sd5_1,sd6_1,sd7_1,sd8_1,sd9_1,sd10_1,sd11_1,sd12_1,sd13_1,sd14_1,sd15_1,sd16_1,sd17_1,sd18_1),2)

d_comparison_1 = data.frame(colnames,mean_d_original,mean_d,sd_1)
d_comparison_1

mean(d_comparison_1$sd_1)

stargazer(d_comparison_1, 
          type = "text",
          colnames = TRUE,
          summary = FALSE,
          title="", 
          digits=2, 
          rownames=FALSE, 
          float = TRUE, 
          float.env = "table", 
          table.placement = "H",
          out="~/Dropbox/Crime Chile/11_replication/tableA3.html")

############################################
# Standardized differences d and d1
############################################

sd1_2 = abs((mean(d$capitalprovince, na.rm = T) - mean(d1$capitalprovince, na.rm = T))/ sqrt( (sd(d$capitalprovince, na.rm = T)^2 - sd(d1$capitalprovince, na.rm = T)^2/2) ))
sd2_2 = abs((mean(d$codeprovince, na.rm = T) - mean(d1$codeprovince, na.rm = T))/ sqrt( (sd(d$codeprovince, na.rm = T)^2 - sd(d1$codeprovince, na.rm = T)^2/2) ))
sd3_2 = abs((mean(d$coderegion, na.rm = T) - mean(d1$coderegion, na.rm = T))/ sqrt( (sd(d$coderegion, na.rm = T)^2 - sd(d1$coderegion, na.rm = T)^2/2) ))
sd4_2 = abs((mean(d$area, na.rm = T) - mean(d1$area, na.rm = T))/ sqrt( (sd(d$area, na.rm = T)^2 - sd(d1$area, na.rm = T)^2/2) ))
sd5_2 = abs((mean(d$logpop02, na.rm = T) - mean(d1$logpop02, na.rm = T))/ sqrt( (sd(d$logpop02, na.rm = T)^2 - sd(d1$logpop02, na.rm = T)^2/2) ))
sd6_2 = abs((mean(d$urbano02_p, na.rm = T) - mean(d1$urbano02_p, na.rm = T))/ sqrt( (sd(d$urbano02_p, na.rm = T)^2 - sd(d1$urbano02_p, na.rm = T)^2/2) ))
sd7_2 = abs((mean(d$alfab03, na.rm = T) - mean(d1$alfab03, na.rm = T))/ sqrt( (sd(d$alfab03, na.rm = T)^2 - sd(d1$alfab03, na.rm = T)^2/2) ))
sd8_2 = abs((mean(d$ingresopercapita03, na.rm = T) - mean(d1$ingresopercapita03, na.rm = T))/ sqrt( (sd(d$ingresopercapita03, na.rm = T)^2 - sd(d1$ingresopercapita03, na.rm = T)^2/2) ))
sd9_2 = abs((mean(d$idh03, na.rm = T) - mean(d1$idh03, na.rm = T))/ sqrt( (sd(d$idh03, na.rm = T)^2 - sd(d1$idh03, na.rm = T)^2/2) ))
sd10_2 = abs((mean(d$rankidh03, na.rm = T) - mean(d1$rankidh03, na.rm = T))/ sqrt( (sd(d$rankidh03, na.rm = T)^2 - sd(d1$rankidh03, na.rm = T)^2/2) ))
sd11_2 = abs((mean(d$logpresvotes1999_1_all, na.rm = T) - mean(d1$logpresvotes1999_1_all, na.rm = T))/ sqrt( (sd(d$logpresvotes1999_1_all, na.rm = T)^2 - sd(d1$logpresvotes1999_1_all, na.rm = T)^2/2) ))
sd12_2 = abs((mean(d$lavin99_p, na.rm = T) - mean(d1$lavin99_p, na.rm = T))/ sqrt( (sd(d$lavin99_p, na.rm = T)^2 - sd(d1$lavin99_p, na.rm = T)^2/2) ))
sd13_2 = abs((mean(d$lagos99_p, na.rm = T) - mean(d1$lagos99_p, na.rm = T))/ sqrt( (sd(d$lagos99_p, na.rm = T)^2 - sd(d1$lagos99_p, na.rm = T)^2/2) ))
sd14_2 = abs((mean(d$nulo99_p, na.rm = T) - mean(d1$nulo99_p, na.rm = T))/ sqrt( (sd(d$nulo99_p, na.rm = T)^2 - sd(d1$nulo99_p, na.rm = T)^2/2) ))
sd15_2 = abs((mean(d$blanco99_p, na.rm = T) - mean(d1$blanco99_p, na.rm = T))/ sqrt( (sd(d$blanco99_p, na.rm = T)^2 - sd(d1$blanco99_p, na.rm = T)^2/2) ))
sd16_2 = abs((mean(d$unemployment_2003, na.rm = T) - mean(d1$unemployment_2003, na.rm = T))/ sqrt( (sd(d$unemployment_2003, na.rm = T)^2 - sd(d1$unemployment_2003, na.rm = T)^2/2) ))
sd17_2 = abs((mean(d$wage_2003, na.rm = T) - mean(d1$wage_2003, na.rm = T))/ sqrt( (sd(d$wage_2003, na.rm = T)^2 - sd(d1$wage_2003, na.rm = T)^2/2) ))
sd18_2 = abs((mean(d$inequality_2003, na.rm = T) - mean(d1$inequality_2003, na.rm = T))/ sqrt( (sd(d$inequality_2003, na.rm = T)^2 - sd(d1$inequality_2003, na.rm = T)^2/2) ))

sd_2 = round(c(sd1_2,sd2_2,sd3_2,sd4_2,sd5_2,sd6_2,sd7_2,sd8_2,sd9_2,sd10_2,sd11_2,sd12_2,sd13_2,sd14_2,sd15_2,sd16_2,sd17_2,sd18_2),2)

d_comparison_2 = data.frame(colnames,mean_d,mean_d1,sd_2)
d_comparison_2

mean(d_comparison_2$sd_2)

stargazer(d_comparison_2, 
          type = "text",
          colnames = TRUE,
          summary = FALSE,
          title="", 
          digits=2, 
          rownames=FALSE, 
          float = TRUE, 
          float.env = "table", 
          table.placement = "H",
          out="~/Dropbox/Crime Chile/11_replication/tableA4.html")

############################################
# Standardized differences d and d2
############################################

sd1_3 = abs((mean(d$capitalprovince, na.rm = T) - mean(d2$capitalprovince, na.rm = T))/ sqrt( (sd(d$capitalprovince, na.rm = T)^2 - sd(d2$capitalprovince, na.rm = T)^2/2) ))
sd2_3 = abs((mean(d$codeprovince, na.rm = T) - mean(d2$codeprovince, na.rm = T))/ sqrt( (sd(d$codeprovince, na.rm = T)^2 - sd(d2$codeprovince, na.rm = T)^2/2) ))
sd3_3 = abs((mean(d$coderegion, na.rm = T) - mean(d2$coderegion, na.rm = T))/ sqrt( (sd(d$coderegion, na.rm = T)^2 - sd(d2$coderegion, na.rm = T)^2/2) ))
sd4_3 = abs((mean(d$area, na.rm = T) - mean(d2$area, na.rm = T))/ sqrt( (sd(d$area, na.rm = T)^2 - sd(d2$area, na.rm = T)^2/2) ))
sd5_3 = abs((mean(d$logpop02, na.rm = T) - mean(d2$logpop02, na.rm = T))/ sqrt( (sd(d$logpop02, na.rm = T)^2 - sd(d2$logpop02, na.rm = T)^2/2) ))
sd6_3 = abs((mean(d$urbano02_p, na.rm = T) - mean(d2$urbano02_p, na.rm = T))/ sqrt( (sd(d$urbano02_p, na.rm = T)^2 - sd(d2$urbano02_p, na.rm = T)^2/2) ))
sd7_3 = abs((mean(d$alfab03, na.rm = T) - mean(d2$alfab03, na.rm = T))/ sqrt( (sd(d$alfab03, na.rm = T)^2 - sd(d2$alfab03, na.rm = T)^2/2) ))
sd8_3 = abs((mean(d$ingresopercapita03, na.rm = T) - mean(d2$ingresopercapita03, na.rm = T))/ sqrt( (sd(d$ingresopercapita03, na.rm = T)^2 - sd(d2$ingresopercapita03, na.rm = T)^2/2) ))
sd9_3 = abs((mean(d$idh03, na.rm = T) - mean(d2$idh03, na.rm = T))/ sqrt( (sd(d$idh03, na.rm = T)^2 - sd(d2$idh03, na.rm = T)^2/2) ))
sd10_3 = abs((mean(d$rankidh03, na.rm = T) - mean(d2$rankidh03, na.rm = T))/ sqrt( (sd(d$rankidh03, na.rm = T)^2 - sd(d2$rankidh03, na.rm = T)^2/2) ))
sd11_3 = abs((mean(d$logpresvotes1999_1_all, na.rm = T) - mean(d2$logpresvotes1999_1_all, na.rm = T))/ sqrt( (sd(d$logpresvotes1999_1_all, na.rm = T)^2 - sd(d2$logpresvotes1999_1_all, na.rm = T)^2/2) ))
sd12_3 = abs((mean(d$lavin99_p, na.rm = T) - mean(d2$lavin99_p, na.rm = T))/ sqrt( (sd(d$lavin99_p, na.rm = T)^2 - sd(d2$lavin99_p, na.rm = T)^2/2) ))
sd13_3 = abs((mean(d$lagos99_p, na.rm = T) - mean(d2$lagos99_p, na.rm = T))/ sqrt( (sd(d$lagos99_p, na.rm = T)^2 - sd(d2$lagos99_p, na.rm = T)^2/2) ))
sd14_3 = abs((mean(d$nulo99_p, na.rm = T) - mean(d2$nulo99_p, na.rm = T))/ sqrt( (sd(d$nulo99_p, na.rm = T)^2 - sd(d2$nulo99_p, na.rm = T)^2/2) ))
sd15_3 = abs((mean(d$blanco99_p, na.rm = T) - mean(d2$blanco99_p, na.rm = T))/ sqrt( (sd(d$blanco99_p, na.rm = T)^2 - sd(d2$blanco99_p, na.rm = T)^2/2) ))
sd16_3 = abs((mean(d$unemployment_2003, na.rm = T) - mean(d2$unemployment_2003, na.rm = T))/ sqrt( (sd(d$unemployment_2003, na.rm = T)^2 - sd(d2$unemployment_2003, na.rm = T)^2/2) ))
sd17_3 = abs((mean(d$wage_2003, na.rm = T) - mean(d2$wage_2003, na.rm = T))/ sqrt( (sd(d$wage_2003, na.rm = T)^2 - sd(d2$wage_2003, na.rm = T)^2/2) ))
sd18_3 = abs((mean(d$inequality_2003, na.rm = T) - mean(d2$inequality_2003, na.rm = T))/ sqrt( (sd(d$inequality_2003, na.rm = T)^2 - sd(d2$inequality_2003, na.rm = T)^2/2) ))

sd_3 = round(c(sd1_3,sd2_3,sd3_3,sd4_3,sd5_3,sd6_3,sd7_3,sd8_3,sd9_3,sd10_3,sd11_3,sd12_3,sd13_3,sd14_3,sd15_3,sd16_3,sd17_3,sd18_3),2)

d_comparison_3 = data.frame(colnames,mean_d,mean_d2,sd_3)
d_comparison_3

mean(d_comparison_3$sd_3)

stargazer(d_comparison_3, 
          type = "text",
          colnames = TRUE,
          summary = FALSE,
          title="", 
          digits=2, 
          rownames=FALSE, 
          float = TRUE, 
          float.env = "table", 
          table.placement = "H",
          out="~/Dropbox/Crime Chile/11_replication/tableA5.html")

############################################
# Standardized differences d and d3
############################################

sd1_4 = abs((mean(d$capitalprovince, na.rm = T) - mean(d3$capitalprovince, na.rm = T))/ sqrt( (sd(d$capitalprovince, na.rm = T)^2 - sd(d3$capitalprovince, na.rm = T)^2/2) ))
sd2_4 = abs((mean(d$codeprovince, na.rm = T) - mean(d3$codeprovince, na.rm = T))/ sqrt( (sd(d$codeprovince, na.rm = T)^2 - sd(d3$codeprovince, na.rm = T)^2/2) ))
sd3_4 = abs((mean(d$coderegion, na.rm = T) - mean(d3$coderegion, na.rm = T))/ sqrt( (sd(d$coderegion, na.rm = T)^2 - sd(d3$coderegion, na.rm = T)^2/2) ))
sd4_4 = abs((mean(d$area, na.rm = T) - mean(d3$area, na.rm = T))/ sqrt( (sd(d$area, na.rm = T)^2 - sd(d3$area, na.rm = T)^2/2) ))
sd5_4 = abs((mean(d$logpop02, na.rm = T) - mean(d3$logpop02, na.rm = T))/ sqrt( (sd(d$logpop02, na.rm = T)^2 - sd(d3$logpop02, na.rm = T)^2/2) ))
sd6_4 = abs((mean(d$urbano02_p, na.rm = T) - mean(d3$urbano02_p, na.rm = T))/ sqrt( (sd(d$urbano02_p, na.rm = T)^2 - sd(d3$urbano02_p, na.rm = T)^2/2) ))
sd7_4 = abs((mean(d$alfab03, na.rm = T) - mean(d3$alfab03, na.rm = T))/ sqrt( (sd(d$alfab03, na.rm = T)^2 - sd(d3$alfab03, na.rm = T)^2/2) ))
sd8_4 = abs((mean(d$ingresopercapita03, na.rm = T) - mean(d3$ingresopercapita03, na.rm = T))/ sqrt( (sd(d$ingresopercapita03, na.rm = T)^2 - sd(d3$ingresopercapita03, na.rm = T)^2/2) ))
sd9_4 = abs((mean(d$idh03, na.rm = T) - mean(d3$idh03, na.rm = T))/ sqrt( (sd(d$idh03, na.rm = T)^2 - sd(d3$idh03, na.rm = T)^2/2) ))
sd10_4 = abs((mean(d$rankidh03, na.rm = T) - mean(d3$rankidh03, na.rm = T))/ sqrt( (sd(d$rankidh03, na.rm = T)^2 - sd(d3$rankidh03, na.rm = T)^2/2) ))
sd11_4 = abs((mean(d$logpresvotes1999_1_all, na.rm = T) - mean(d3$logpresvotes1999_1_all, na.rm = T))/ sqrt( (sd(d$logpresvotes1999_1_all, na.rm = T)^2 - sd(d3$logpresvotes1999_1_all, na.rm = T)^2/2) ))
sd12_4 = abs((mean(d$lavin99_p, na.rm = T) - mean(d3$lavin99_p, na.rm = T))/ sqrt( (sd(d$lavin99_p, na.rm = T)^2 - sd(d3$lavin99_p, na.rm = T)^2/2) ))
sd13_4 = abs((mean(d$lagos99_p, na.rm = T) - mean(d3$lagos99_p, na.rm = T))/ sqrt( (sd(d$lagos99_p, na.rm = T)^2 - sd(d3$lagos99_p, na.rm = T)^2/2) ))
sd14_4 = abs((mean(d$nulo99_p, na.rm = T) - mean(d3$nulo99_p, na.rm = T))/ sqrt( (sd(d$nulo99_p, na.rm = T)^2 - sd(d3$nulo99_p, na.rm = T)^2/2) ))
sd15_4 = abs((mean(d$blanco99_p, na.rm = T) - mean(d3$blanco99_p, na.rm = T))/ sqrt( (sd(d$blanco99_p, na.rm = T)^2 - sd(d3$blanco99_p, na.rm = T)^2/2) ))
sd16_4 = abs((mean(d$unemployment_2003, na.rm = T) - mean(d3$unemployment_2003, na.rm = T))/ sqrt( (sd(d$unemployment_2003, na.rm = T)^2 - sd(d3$unemployment_2003, na.rm = T)^2/2) ))
sd17_4 = abs((mean(d$wage_2003, na.rm = T) - mean(d3$wage_2003, na.rm = T))/ sqrt( (sd(d$wage_2003, na.rm = T)^2 - sd(d3$wage_2003, na.rm = T)^2/2) ))
sd18_4 = abs((mean(d$inequality_2003, na.rm = T) - mean(d3$inequality_2003, na.rm = T))/ sqrt( (sd(d$inequality_2003, na.rm = T)^2 - sd(d3$inequality_2003, na.rm = T)^2/2) ))

sd_4 = round(c(sd1_4,sd2_4,sd3_4,sd4_4,sd5_4,sd6_4,sd7_4,sd8_4,sd9_4,sd10_4,sd11_4,sd12_4,sd13_4,sd14_4,sd15_4,sd16_4,sd17_4,sd18_4),2)

d_comparison_4 = data.frame(colnames,mean_d,mean_d3,sd_4)
d_comparison_4

mean(d_comparison_4$sd_4)

stargazer(d_comparison_4, 
          type = "text",
          colnames = TRUE,
          summary = FALSE,
          title="", 
          digits=2, 
          rownames=FALSE, 
          float = TRUE, 
          float.env = "table", 
          table.placement = "H",
          out="~/Dropbox/Crime Chile/11_replication/tableA6.html")

############################################
# Standardized differences d and d4
############################################

sd1_5 = abs((mean(d$capitalprovince, na.rm = T) - mean(d4$capitalprovince, na.rm = T))/ sqrt( (sd(d$capitalprovince, na.rm = T)^2 - sd(d4$capitalprovince, na.rm = T)^2/2) ))
sd2_5 = abs((mean(d$codeprovince, na.rm = T) - mean(d4$codeprovince, na.rm = T))/ sqrt( (sd(d$codeprovince, na.rm = T)^2 - sd(d4$codeprovince, na.rm = T)^2/2) ))
sd3_5 = abs((mean(d$coderegion, na.rm = T) - mean(d4$coderegion, na.rm = T))/ sqrt( (sd(d$coderegion, na.rm = T)^2 - sd(d4$coderegion, na.rm = T)^2/2) ))
sd4_5 = abs((mean(d$area, na.rm = T) - mean(d4$area, na.rm = T))/ sqrt( (sd(d$area, na.rm = T)^2 - sd(d4$area, na.rm = T)^2/2) ))
sd5_5 = abs((mean(d$logpop02, na.rm = T) - mean(d4$logpop02, na.rm = T))/ sqrt( (sd(d$logpop02, na.rm = T)^2 - sd(d4$logpop02, na.rm = T)^2/2) ))
sd6_5 = abs((mean(d$urbano02_p, na.rm = T) - mean(d4$urbano02_p, na.rm = T))/ sqrt( (sd(d$urbano02_p, na.rm = T)^2 - sd(d4$urbano02_p, na.rm = T)^2/2) ))
sd7_5 = abs((mean(d$alfab03, na.rm = T) - mean(d4$alfab03, na.rm = T))/ sqrt( (sd(d$alfab03, na.rm = T)^2 - sd(d4$alfab03, na.rm = T)^2/2) ))
sd8_5 = abs((mean(d$ingresopercapita03, na.rm = T) - mean(d4$ingresopercapita03, na.rm = T))/ sqrt( (sd(d$ingresopercapita03, na.rm = T)^2 - sd(d4$ingresopercapita03, na.rm = T)^2/2) ))
sd9_5 = abs((mean(d$idh03, na.rm = T) - mean(d4$idh03, na.rm = T))/ sqrt( (sd(d$idh03, na.rm = T)^2 - sd(d4$idh03, na.rm = T)^2/2) ))
sd10_5 = abs((mean(d$rankidh03, na.rm = T) - mean(d4$rankidh03, na.rm = T))/ sqrt( (sd(d$rankidh03, na.rm = T)^2 - sd(d4$rankidh03, na.rm = T)^2/2) ))
sd11_5 = abs((mean(d$logpresvotes1999_1_all, na.rm = T) - mean(d4$logpresvotes1999_1_all, na.rm = T))/ sqrt( (sd(d$logpresvotes1999_1_all, na.rm = T)^2 - sd(d4$logpresvotes1999_1_all, na.rm = T)^2/2) ))
sd12_5 = abs((mean(d$lavin99_p, na.rm = T) - mean(d4$lavin99_p, na.rm = T))/ sqrt( (sd(d$lavin99_p, na.rm = T)^2 - sd(d4$lavin99_p, na.rm = T)^2/2) ))
sd13_5 = abs((mean(d$lagos99_p, na.rm = T) - mean(d4$lagos99_p, na.rm = T))/ sqrt( (sd(d$lagos99_p, na.rm = T)^2 - sd(d4$lagos99_p, na.rm = T)^2/2) ))
sd14_5 = abs((mean(d$nulo99_p, na.rm = T) - mean(d4$nulo99_p, na.rm = T))/ sqrt( (sd(d$nulo99_p, na.rm = T)^2 - sd(d4$nulo99_p, na.rm = T)^2/2) ))
sd15_5 = abs((mean(d$blanco99_p, na.rm = T) - mean(d4$blanco99_p, na.rm = T))/ sqrt( (sd(d$blanco99_p, na.rm = T)^2 - sd(d4$blanco99_p, na.rm = T)^2/2) ))
sd16_5 = abs((mean(d$unemployment_2003, na.rm = T) - mean(d4$unemployment_2003, na.rm = T))/ sqrt( (sd(d$unemployment_2003, na.rm = T)^2 - sd(d4$unemployment_2003, na.rm = T)^2/2) ))
sd17_5 = abs((mean(d$wage_2003, na.rm = T) - mean(d4$wage_2003, na.rm = T))/ sqrt( (sd(d$wage_2003, na.rm = T)^2 - sd(d4$wage_2003, na.rm = T)^2/2) ))
sd18_5 = abs((mean(d$inequality_2003, na.rm = T) - mean(d4$inequality_2003, na.rm = T))/ sqrt( (sd(d$inequality_2003, na.rm = T)^2 - sd(d4$inequality_2003, na.rm = T)^2/2) ))

sd_5 = round(c(sd1_5,sd2_5,sd3_5,sd4_5,sd5_5,sd6_5,sd7_5,sd8_5,sd9_5,sd10_5,sd11_5,sd12_5,sd13_5,sd14_5,sd15_5,sd16_5,sd17_5,sd18_5),2)

d_comparison_5 = data.frame(colnames,mean_d,mean_d4,sd_5)
d_comparison_5

mean(d_comparison_5$sd_5)

stargazer(d_comparison_5, 
          type = "text",
          colnames = TRUE,
          summary = FALSE,
          title="", 
          digits=2, 
          rownames=FALSE, 
          float = TRUE, 
          float.env = "table", 
          table.placement = "H",
          out="~/Dropbox/Crime Chile/11_replication/tableA7.html")

############################################
# Standardized differences d and d5
############################################

sd1_6 = abs((mean(d$capitalprovince, na.rm = T) - mean(d5$capitalprovince, na.rm = T))/ sqrt( (sd(d$capitalprovince, na.rm = T)^2 - sd(d5$capitalprovince, na.rm = T)^2/2) ))
sd2_6 = abs((mean(d$codeprovince, na.rm = T) - mean(d5$codeprovince, na.rm = T))/ sqrt( (sd(d$codeprovince, na.rm = T)^2 - sd(d5$codeprovince, na.rm = T)^2/2) ))
sd3_6 = abs((mean(d$coderegion, na.rm = T) - mean(d5$coderegion, na.rm = T))/ sqrt( (sd(d$coderegion, na.rm = T)^2 - sd(d5$coderegion, na.rm = T)^2/2) ))
sd4_6 = abs((mean(d$area, na.rm = T) - mean(d5$area, na.rm = T))/ sqrt( (sd(d$area, na.rm = T)^2 - sd(d5$area, na.rm = T)^2/2) ))
sd5_6 = abs((mean(d$logpop02, na.rm = T) - mean(d5$logpop02, na.rm = T))/ sqrt( (sd(d$logpop02, na.rm = T)^2 - sd(d5$logpop02, na.rm = T)^2/2) ))
sd6_6 = abs((mean(d$urbano02_p, na.rm = T) - mean(d5$urbano02_p, na.rm = T))/ sqrt( (sd(d$urbano02_p, na.rm = T)^2 - sd(d5$urbano02_p, na.rm = T)^2/2) ))
sd7_6 = abs((mean(d$alfab03, na.rm = T) - mean(d5$alfab03, na.rm = T))/ sqrt( (sd(d$alfab03, na.rm = T)^2 - sd(d5$alfab03, na.rm = T)^2/2) ))
sd8_6 = abs((mean(d$ingresopercapita03, na.rm = T) - mean(d5$ingresopercapita03, na.rm = T))/ sqrt( (sd(d$ingresopercapita03, na.rm = T)^2 - sd(d5$ingresopercapita03, na.rm = T)^2/2) ))
sd9_6 = abs((mean(d$idh03, na.rm = T) - mean(d5$idh03, na.rm = T))/ sqrt( (sd(d$idh03, na.rm = T)^2 - sd(d5$idh03, na.rm = T)^2/2) ))
sd10_6 = abs((mean(d$rankidh03, na.rm = T) - mean(d5$rankidh03, na.rm = T))/ sqrt( (sd(d$rankidh03, na.rm = T)^2 - sd(d5$rankidh03, na.rm = T)^2/2) ))
sd11_6 = abs((mean(d$logpresvotes1999_1_all, na.rm = T) - mean(d5$logpresvotes1999_1_all, na.rm = T))/ sqrt( (sd(d$logpresvotes1999_1_all, na.rm = T)^2 - sd(d5$logpresvotes1999_1_all, na.rm = T)^2/2) ))
sd12_6 = abs((mean(d$lavin99_p, na.rm = T) - mean(d5$lavin99_p, na.rm = T))/ sqrt( (sd(d$lavin99_p, na.rm = T)^2 - sd(d5$lavin99_p, na.rm = T)^2/2) ))
sd13_6 = abs((mean(d$lagos99_p, na.rm = T) - mean(d5$lagos99_p, na.rm = T))/ sqrt( (sd(d$lagos99_p, na.rm = T)^2 - sd(d5$lagos99_p, na.rm = T)^2/2) ))
sd14_6 = abs((mean(d$nulo99_p, na.rm = T) - mean(d5$nulo99_p, na.rm = T))/ sqrt( (sd(d$nulo99_p, na.rm = T)^2 - sd(d5$nulo99_p, na.rm = T)^2/2) ))
sd15_6 = abs((mean(d$blanco99_p, na.rm = T) - mean(d5$blanco99_p, na.rm = T))/ sqrt( (sd(d$blanco99_p, na.rm = T)^2 - sd(d5$blanco99_p, na.rm = T)^2/2) ))
sd16_6 = abs((mean(d$unemployment_2003, na.rm = T) - mean(d5$unemployment_2003, na.rm = T))/ sqrt( (sd(d$unemployment_2003, na.rm = T)^2 - sd(d5$unemployment_2003, na.rm = T)^2/2) ))
sd17_6 = abs((mean(d$wage_2003, na.rm = T) - mean(d5$wage_2003, na.rm = T))/ sqrt( (sd(d$wage_2003, na.rm = T)^2 - sd(d5$wage_2003, na.rm = T)^2/2) ))
sd18_6 = abs((mean(d$inequality_2003, na.rm = T) - mean(d5$inequality_2003, na.rm = T))/ sqrt( (sd(d$inequality_2003, na.rm = T)^2 - sd(d5$inequality_2003, na.rm = T)^2/2) ))

sd_6 = round(c(sd1_6,sd2_6,sd3_6,sd4_6,sd5_6,sd6_6,sd7_6,sd8_6,sd9_6,sd10_6,sd11_6,sd12_6,sd13_6,sd14_6,sd15_6,sd16_6,sd17_6,sd18_6),2)

d_comparison_6 = data.frame(colnames,mean_d,mean_d5,sd_6)
d_comparison_6

mean(d_comparison_6$sd_6)

stargazer(d_comparison_6, 
          type = "text",
          colnames = TRUE,
          summary = FALSE,
          title="", 
          digits=2, 
          rownames=FALSE, 
          float = TRUE, 
          float.env = "table", 
          table.placement = "H",
          out="~/Dropbox/Crime Chile/11_replication/tableA8.html")

############################################
# Standardized differences d and d5
############################################

sd1_7 = abs((mean(d$capitalprovince, na.rm = T) - mean(d6$capitalprovince, na.rm = T))/ sqrt( (sd(d$capitalprovince, na.rm = T)^2 - sd(d6$capitalprovince, na.rm = T)^2/2) ))
sd2_7 = abs((mean(d$codeprovince, na.rm = T) - mean(d6$codeprovince, na.rm = T))/ sqrt( (sd(d$codeprovince, na.rm = T)^2 - sd(d6$codeprovince, na.rm = T)^2/2) ))
sd3_7 = abs((mean(d$coderegion, na.rm = T) - mean(d6$coderegion, na.rm = T))/ sqrt( (sd(d$coderegion, na.rm = T)^2 - sd(d6$coderegion, na.rm = T)^2/2) ))
sd4_7 = abs((mean(d$area, na.rm = T) - mean(d6$area, na.rm = T))/ sqrt( (sd(d$area, na.rm = T)^2 - sd(d6$area, na.rm = T)^2/2) ))
sd5_7 = abs((mean(d$logpop02, na.rm = T) - mean(d6$logpop02, na.rm = T))/ sqrt( (sd(d$logpop02, na.rm = T)^2 - sd(d6$logpop02, na.rm = T)^2/2) ))
sd6_7 = abs((mean(d$urbano02_p, na.rm = T) - mean(d6$urbano02_p, na.rm = T))/ sqrt( (sd(d$urbano02_p, na.rm = T)^2 - sd(d6$urbano02_p, na.rm = T)^2/2) ))
sd7_7 = abs((mean(d$alfab03, na.rm = T) - mean(d6$alfab03, na.rm = T))/ sqrt( (sd(d$alfab03, na.rm = T)^2 - sd(d6$alfab03, na.rm = T)^2/2) ))
sd8_7 = abs((mean(d$ingresopercapita03, na.rm = T) - mean(d6$ingresopercapita03, na.rm = T))/ sqrt( (sd(d$ingresopercapita03, na.rm = T)^2 - sd(d6$ingresopercapita03, na.rm = T)^2/2) ))
sd9_7 = abs((mean(d$idh03, na.rm = T) - mean(d6$idh03, na.rm = T))/ sqrt( (sd(d$idh03, na.rm = T)^2 - sd(d6$idh03, na.rm = T)^2/2) ))
sd10_7 = abs((mean(d$rankidh03, na.rm = T) - mean(d6$rankidh03, na.rm = T))/ sqrt( (sd(d$rankidh03, na.rm = T)^2 - sd(d6$rankidh03, na.rm = T)^2/2) ))
sd11_7 = abs((mean(d$logpresvotes1999_1_all, na.rm = T) - mean(d6$logpresvotes1999_1_all, na.rm = T))/ sqrt( (sd(d$logpresvotes1999_1_all, na.rm = T)^2 - sd(d6$logpresvotes1999_1_all, na.rm = T)^2/2) ))
sd12_7 = abs((mean(d$lavin99_p, na.rm = T) - mean(d6$lavin99_p, na.rm = T))/ sqrt( (sd(d$lavin99_p, na.rm = T)^2 - sd(d6$lavin99_p, na.rm = T)^2/2) ))
sd13_7 = abs((mean(d$lagos99_p, na.rm = T) - mean(d6$lagos99_p, na.rm = T))/ sqrt( (sd(d$lagos99_p, na.rm = T)^2 - sd(d6$lagos99_p, na.rm = T)^2/2) ))
sd14_7 = abs((mean(d$nulo99_p, na.rm = T) - mean(d6$nulo99_p, na.rm = T))/ sqrt( (sd(d$nulo99_p, na.rm = T)^2 - sd(d6$nulo99_p, na.rm = T)^2/2) ))
sd15_7 = abs((mean(d$blanco99_p, na.rm = T) - mean(d6$blanco99_p, na.rm = T))/ sqrt( (sd(d$blanco99_p, na.rm = T)^2 - sd(d6$blanco99_p, na.rm = T)^2/2) ))
sd16_7 = abs((mean(d$unemployment_2003, na.rm = T) - mean(d6$unemployment_2003, na.rm = T))/ sqrt( (sd(d$unemployment_2003, na.rm = T)^2 - sd(d6$unemployment_2003, na.rm = T)^2/2) ))
sd17_7 = abs((mean(d$wage_2003, na.rm = T) - mean(d6$wage_2003, na.rm = T))/ sqrt( (sd(d$wage_2003, na.rm = T)^2 - sd(d6$wage_2003, na.rm = T)^2/2) ))
sd18_7 = abs((mean(d$inequality_2003, na.rm = T) - mean(d6$inequality_2003, na.rm = T))/ sqrt( (sd(d$inequality_2003, na.rm = T)^2 - sd(d6$inequality_2003, na.rm = T)^2/2) ))

sd_7 = round(c(sd1_7,sd2_7,sd3_7,sd4_7,sd5_7,sd6_7,sd7_7,sd8_7,sd9_7,sd10_7,sd11_7,sd12_7,sd13_7,sd14_7,sd15_7,sd16_7,sd17_7,sd18_7),2)

d_comparison_7 = data.frame(colnames,mean_d,mean_d6,sd_7)
d_comparison_7

mean(d_comparison_7$sd_7)

stargazer(d_comparison_7, 
          type = "text",
          colnames = TRUE,
          summary = FALSE,
          title="", 
          digits=2, 
          rownames=FALSE, 
          float = TRUE, 
          float.env = "table", 
          table.placement = "H",
          out="~/Dropbox/Crime Chile/11_replication/tableA9.html")

sink()
